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Abstract. We consider the discontinuous Galerkin 

laws derivable as moments of the kmetjc ^ n J linear conser vation laws equipped with an entropy extension 
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1 Introduction Discontinuous Galerkin (DG) finite element methods for first order hyperbolic equa- 

law fluxes ambiguous thus necessitating the introduction of a numerical flux function, ( - > + - {* , d 

yS p Barth [2 ll This latter energy technique is used in the present analysis. Using the notation intro uce 
SS S, initial- value problem one obtains fro. 
energy balance equation for the DG finite element method for a spatial domain $7 integrated 
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This exact energy balance is derived using either of two different baseline numerical flux functions: 

• Symmetric Mean- Value (SMV) Flux 


b S Mv(v-,v + ; n) = i(f(v_) + f(v + )) - ~ hsMv( v -> v +; n ) 


with 


bsMv( v -> v +;n) = jT l |A(vW;n)U o ^[v]t . 


• Kinetic Symmetric Mean- Value (KSMV) Flux 


h K SMv(v-,v + ;n) = i(f(v_) + f(v + )) - ^hKSMv( v -. v +;n) 
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with 

z * 1 

hKSMv( v -. v +; n ) = / (|vn|m®mexp(v(0) • m(u))><W[v]I 
Jo 

where (•) denotes an integration in velocity-internal energy phase space and m(u) the vector of moments as 
discussed later. Observe that the nonlinear energy balance (1.1) formally bounds the final solution in terms 
of initial data. The discontinuous function space leads to energy removal in space and time proportional to 
the matrix modulated square of solution jumps across the respective space and time interfaces. 

In general, the numerical flux functions given here are too complex to permit the calculation of the needed 
path integrations in closed form. Our strategy in this paper is to first develop the framework and theoretical 
results given above. We then shift to our main objective, the development of approximate numerical flux 
functions, h ap prox(v-, v+; n), that avoid these complicated path integrations without compromising our 
ability to rigorously prove nonlinear stability. This is accomplished by requiring that the approximate 
numerical flux formulas are more energy dissipative than the theoretically derived fluxes given above. is 
task can be reduced to the satisfaction of either of two algebraic sufficient conditions (derived later) 


[V]£ • hg MV < [*]£ ' hfpprox « r Ml- ' h KSMV < ML ' h 


d 

approx 


We then construct a number of approximate flux functions based on this strategy. 

2. Background. Consider the Cauchy initial value problem for a system of m coupled first-order 
differential equations in d space coordinates and time which represents a conservation law process. Let 
u(x,t) : IR d x 1R + t-4 lR m denote the dependent solution variables and f(u) : lR m *-4 IR mx the flux vector. 
The prototype Cauchy problem is then given by 

/ u t + r Xi = o 

(2*1) \ u(x, 0) = uo(^) 

with implied summation on the index t. Additionally, the system is assumed to possess an scalar entropy 
extension. Let U{ u) : IR m i-4 IR denote an entropy function and F(u) : H m »-> 1R the entropy flux such 
that in addition to (2.1) the following inequality holds 


(2.2) U,t + F 'x< - 0 

with equality for smooth solutions. In symmetrization theory for first-order conservation laws [11, 17, 12] , 
one seeks a mapping u(v) : IR m >-4 IR m applied to (2.1) so that when transformed 

(2.3) u ,vV,t + r v v iXi = 0 

the matrix u, v is symmetric positive definite (SPD) and the matrices T v are symmetric. Clearly, if functions 
U(y) : H m t-4 IR and F'(v) : lR m t-4 IR can be found so that 

(2.4) u T = W.v, (C) T = 


then the matrices 

(2.5) ^. v “ W,v,v> C v *^? v , v 

are symmetric. Further, we shall require that W(v) be a differentiable convex function such that 


( 2 . 6 ) 


lim 

v— >oo 


W(v) 


= +00 


so that U{ u) can be interpreted as a Legendre transform of U(v) 

{/(u) = sup (v • u - W(v)} . 

V 


(2.7) 
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From (2.6), it follows that 3 v* <E lR m such that v u - W(v) achieves a maximum at v* 

(2 .8) £/(u)=v*-u — W(v') • 

At this maximum u = W, v (V) which can be locally inverted to the form v* = v(u). Elimination of v* in 

(2.8) yields the simplified duality relationship 

( 2 . 9 ) U (u) = v(u) • u — W(v(u)) . 


Differentiation of this expression 


(2.10) t/ u = v ^~ ■f' uTv i» — — v 

gives an explicit formula for the entropy variables v in terms of derivatives of the entropy function C/(u) 

(2.11) vT = ■ 

Using the mapping relation v(u), a duality pairing for entropy flux components is defined 


( 2 . 12 ) 


F’(u) = v(u) • r(u) - ^(v(u)) . 


Differentiation then yields the flux relation 

(2.i3) f!u = v T ru + (r) T v, u 


F>,, 


= v T f*„ 


and the fundamental relationship for smooth solutions 

(2.14) v ' (u,t + f^ Zi ) = U,t + F' x . = 0 


which is exploited in nonlinear energy analysis. . 

Note that convexity of U(v) implies positive definiteness of u v and hyperbolicity of (2.1) [8, 17J, viz., 
that the linear combination f, u (n) = n, f' u has real eigenvalues and a complete set of real-valued eigenvectors 

for all nonzero n 6 TR d . This result follows immediately from the identity 

(u, v )- 1/2 f, U (n)(u, v ) 1/2 = (u v ) - 1 / 2 f,v(n)(u, v )" 1/2 

v V “* f 

symm 


which shows that f, u (n) is similar to a symmetric matrix. 

2.1. Kinetic Boltzmann Entropies. Consider the particular case of moment systems derived from 
the kinetic Boltzmann equation with Levermore’s closure [16]. Boltzmann’s equation is given by 

(2.15) f{x,v,t)'t + v ■ V x f(x,v,t) = C(f)(x,v,t) , 

with f(x,v,t) a nonnegative density function, v € JR d the velocity, and C(f) : IR 1R the collision operator. 
Moment systems are obtained by integrating in velocity space the Boltzmann equation over a vector m(t>) 
of linearly independant polynomials in velocity, 

(2.16) (m/), t + = (mC(/)) , 

where (t/>) denotes the integral of a measurable function rjj over velocity space. Without further assumption, 
the fluxes (tij m /) cannot be expressed as functions of u = (m /) . The closure of the system is performed y 
assuming that the distribution function / has a prescribed form / B = /b(u) given by the minimum entropy 
principle 

(2.17) H[f B ] = min{tf [3] |(jm) = u} , 

where H[g] = (g In g) is Boltzmann’s celebrated H-function. Since H is a convex function the minimization 
problem (2.17) is formally equivalent to 


(2.18) 


/ B = exp(v • m) , 
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where v = v(u) serves as the Lagrange multiplier associated with the constraint {g m) — u or equivalently 
under the closure assumption 

(2.19) u = (m exp(v(u) ■ m)) . 

The moment system (2.16) can now be rewritten as 

(2.20) u ,t + ?*< = r ( u > - 

where 

(2.21) r = (Vi m exp(v(u) • m)) . 

Observe that using the kinetic Boltzmann structure, we have that 

(2.22) ^(v) = (/> = ( ex P( v • m )> 

is a suitable conjugate entropy function and that 

(2.23) U{ u) = ((v(u) • m) exp(v(u) • m) - exp(v(u) • m)) 

is the corresponding entropy function so that the duality relationship (2.9) holds. 

The simplest example of a moment system is obtained by taking m(v) = (l,u, \v\ 2 /2) T corresponding 
to mass, momentum, and kinetic energy. In this instance, the collision integral vanishes identically, r = 0, 
and (2.20) is the well-known system of Euler equations (5 moments) for a monotonic gas. More complex 
systems with 10, 14 or 35 moments have been considered in the literature [16]. In Appendix A, we give the 
corresponding Euler equations moment model for 7-law (polytropic) gases that is achieved by increasing the 
dimension of the phase space to include internal energy I and utilizing the moments m(v, I) = (1, v, |v| /2 + 
jf)T f or (5 _ (1/(7 - 1) — d/2) -1 . In the case of the 7-law gas, one obtains a conjugate entropy function in 

JR" 1 of the form 

(2.24) W(v) = (c(7,d) exp(v • m)) , c(7,d) > 0 

which is still compatible with the desired exponential structure. For brevity, we will omit constants such as 
c(7,d) in our exponential form so that (2.22) may be regarded as an abstract form for (2.24) with suitably 
chosen phase space. From (2.22) it is clear that 

(2.25) ZY.v.v = (m ® m exp(v • m)) 
is SPD, i.e. the following double contraction to a scalar is positive 

(2.26) U , Vi , v> ZiZ j = ((m ■ z) 2 exp(v • m)) > 0 , |z| / 0 . 

Furthermore, (7 UU = W',’ v is also SPD, hence U is also a convex function of u. Consequently every system 
with the considered structure is hyperbolic symmetrizable and has a convex entropy U which is locally 
dissipated. This technique provides one of the simplest proofs of convexity for entropy functions associated 
with first order nonlinear conservation law systems derivable as moment closures of kinetic Boltzmann-like 
equations. In the case of the Euler equations of gasdynamics, the reader should compare this technique with 
the somewhat tedious proofs of convexity given in Refs. [12], [13], [9]. Finally, we mention the following 
general result for kinetic Boltzmann moment hierarchical systems which is used in later development. 

Lemma 2.1. Generalized Convexity of Boltzmann Moment Conjugate Entropies. Let IN — 
{0,1,2,...} denote the set of nonnegative integers. All 2k derivatives of the kinetic Boltzmann moment 
conjugate entropy (2.22) 


U(v) = (exp(v • m)) 
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are SPD for k £ IN 
(2.27) 


d 2k U , 


2k times 


M / o . 


Proof. Successive differentiation of (2.22) 2k times yields the symmetric rank-2fc x m tensor 


^2,28) q 2 k ( v ) = ( m 0 m g> ... 0 m exp(v • rn)) , 

2k times 

followed by contraction to a scalar by a nonzero vector z 6 IR m 

(2.29) |^[z 1 z^z) = ((m • z) 2k exp(v • m)) . 

2k times 

The moment vector m contains m linearly independent polynomials spanning lR m . The condition m(u) z = 0 
for fixed nonzero z and variable v would violate the assumption of linear independence, thus we conclude 
that m(i>) z ^ 0 a.e., namely, except at points of measure zero in the phase space Lebesgue integration. e 
term exp(v • m) is also positive for finite argument values in the phase space integration, hence we conclude 
for nonnegative powers 2k 

(2.30) (( m ' z ) 2k ex P( v ’ m )) > 0 


and the stated lemma. ■ 

2.2. The Eigenvector Scaling Theorem and Generalized Matrix Functions with Respect to 
the A 0 Inner Product. Next, we consider an important algebraic property of right symmetnzable systems 
which is used later in the implementation of the DG scheme. Simplifying upon the previous notation, let 
Aq = u v , A, = P u , Ai = AiA 0 and rewrite (2.3) 


(2.31) 


i 0 v, t + Ai\ tXi = 0 . 


The following theorem states a property of the symmetric matrix A { symmetrized via the symmetric positive 


definite matrix -4 0 - 

Theorem 2.2 (Eigenvector Scaling). 
set of all right symmetrizers: 


Let A € JR nx " be an arbitrary diagonalizable matrix and S the 


S = {B € IR nxn | B SPD, AB symmetric}. 


Further, let R 6 IR nxn denote the right eigenvector matrix which diagonalizes A 

A = RAR~ l 


with r distinct eigenvalues, A = Diag(Ai/ mi x m , , A2/ mj xmji • • • i Wm,xm r )- Then for each B € S there 
exists a symmetric block diagonal matnx T = Diag (T m , xmi , T maXmi , . . . ,r m „ xmr ) that block scales columns 
of R, R — RT, such that 

B = RR t , A = RAR~ l 


which imply that 


AB = R\R t . 


Proof. Omitted, see [2]. ■ 
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This last formula states a congruence relationship since R is not generally orthonormal and l A does not 
represent the eigenvalues of AB. We shall refer to R as containing “entropy scaled eigenvectors^ Note that 
we can consider scalar combinations of A, with the same scaling properties for arbitrary n € IK. , l.e. 


(2.32) 


A(n) = n iAi = R( n) A(n) R T ( n) , A 0 = R( n) fl T (n) 


Wavespeeds associated with the system (2.31) and the direction vector n are given by critical values of 
the Rayleigh quotient 


(2.33) 


(7M n) j _ £ T fi(n) A(n) R r (n) { _ rfA{n) r? e H m , 7 ? = -R£, |£| 7 ^ °> , 

e A 0 £ t-R(n)R T (n)Z VV 


which are simply elements of A(n). For use in later developments, it is useful to define a matrix func- 
tion fi (A) with respect to the Riemannian matrix A 0 with critical values of the Rayleigh quotient given 
by /(Aii),i = 1, . . . ,m. This matrix function takes a particularly simple form as given by the following 

proposition: 

Proposition 2.3. Barth (2, If. Let A 0 denote the SPD right symmetrizer of A such that A = AAo, 
Ao = RR t , and A = RAR~ l - The generalized matrix function f^ 0 {A) ts symmetric and defined canonxca y 

in terms of entropy scaled eigenvectors as 

(2.34) hM) = fl/(A)R r 

where /(A) ts performed componentwise. 

Proof. Assume the desired critical values /(A) and the Rayleigh quotient producing them 


(2.35) 


T7 T /(A)r? ( r Rf{A)R T Z £ T /.4 0 (AK ^jjelR"*, n = R£, |£| #0 


T) ■ T) 


£ T R R T £ 


e r A 0 z 


see also [2, 1]. ■ 

In later sections, the generalized matrix absolute value function | A\ Ao will be required. Using Proposition 
(2.3) stated above 

(2.36) \*\A. ~ -R|A|fi T • 

Finally, observe that using these scaled eigenvectors, R, we have the following equivalent representations 1 
of A and A 0 that are used in later developments: 

m m 

(2.37) A = ^ Aj fj ® fj , A 0 = y; fj ® fj , 

1=1 *=1 


and 

(2.38) 


•=1 


where fi denotes the i-th column of R . 

iThese representations should not be confused with the spectral decomposition of a matrix by orthonormal transform. 
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3. DG Finite Element Method. Let ffc denote a spatial domain composed of nonoverlapping elements 
T{, fi = U Ti, Ti C\Tj — i ^ j and I n =]t n ,£ n+1 [ the n-th time interval. It is useful to also define the 
element set T = {T u T 2 , . . . , T m } and edge set £ = {ei, e 2 , . . . , e\ £ \}. To simplify the exposition, consider a 
single variational formulation with weakly enforced boundary conditions. In the DG formulations (see [15, 3] 
and references therein), functions are discontinuous in space and time, i.e. 

v h = {v h I vf Tx/n € (v k (Tx /”))"*} . 

For ease of exposition, we consider a spatial domain fi which is either periodic in all space dimensions or 
nonperiodic with compactly supported initial data. Consider the first order Cauchy system 

/ u * + f ,x. =0 infi 

‘ \ u(x, 0) = u 0 (x) 

with A(n) = n » Ai and A(n) = n» A { . The DG scheme with weakly imposed boundary conditions in time is 
defined by the following stabilized variational formulation: 

Find v h G V h such that for all v/ h G V h 

(3.2) B(v\w h )GAL=0 


where 


S(v,w)gal = [ [ (-u(v) • w if - f’(v) • w <Xi )dxdt 

Ji^Jq 

+ [ (w (t n _ +l ) ■ u(v(f" +1 )) - w (t n + ) • u(v(£))) dx 

Ju 

+ / V /(w(x_) - w(x+)) •h(v(x.-),v(x + );n)dxdt 
Jln 'G£ Je 


where h denotes a numerical flux function. Throughout, we consider numerical fluxes of the form 
(3.3) h(v_,v + ;n) = ^ (f(v_;n) +f(v + ;n)) - |h d (v_,v + ;n) . 

These fluxes are consistent with the true flux in the sense that f(v;n) = h(v,v;n). 

3,1. DG Nonlinear Energy Analysis. Before presenting the nonlinear energy result, we recall some 
supporting corollaries concerning entropy function/flux jump identities at space-time slab interfaces. Note 
that throughout this section, we utilize the state-space parameterization 

V(0)=v(x_)+0[v]*i 


(similarly across time slab interfaces) for use in state-space path integrations and the interface averaging 
operator 


«V» 


x + — 
X _ — 


v(s-) + v(s+) 

2 


Lemma 3.1. Interface Jump Identities. Barth [2, 1] Let Z(u),Z(v) : ]R m IR be twice differentiable 
functions of their argument satisfying the duality relationship 


(3.4) Z(u) +Z(v) = Z v v . 

The folloxving jump identities hold across interfaces 

(3.5) [Z)%t ~ [ZAlt v(x+) + f\ 1 - 0)[v)tf 2,v.v(v(d)) [v]£ dO = 0 

JO 

m*_ - [Z,v]xt v(x_) - f 6 [v]£- 2, V,v(v(6>)) [v]*+ d6 = 0 . 

Jo 


(3.6) 
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Proof. Omitted, see [2, 1]. ■ 

Corollary 3.2. Temporal Space-Time Slab Interface Identity. Barth [2, 1]. Let t± denote 
a temporal space-time slab interface. The following entropy function jump identity holds across time slab 
interfaces 

(3.7) jf ([U) l t + _ - v t (4) H:) dx + 1||| [vfc \\\% n = 0 
where 

(3.8) III [y]\: |||^ >n = JJ^ 2 (1 -9) [y)l: ■ Mm) Ml! dBdx> O . 

Corollary 3.3. Spatial Space-Time Slab Interface Identity. Barth [2, If. Let x± denote a 
spatial element interface. The following entropy flux jump identity holds across spatial element interfaces 

(3.9) [FXl - ((v T y x : [r]*+ + lj\i-20) Ml: ■ &<?&)) M 1: d6 = 0 . 


Note that in actual numerical calculations, it is desirable to use the variational form given by (3.2) since 
integration by parts has been used to insure exact discrete conservation even with inexact numerical quadra- 
ture of the various integrals. For analysis purposes, however, it is desirable to use the following equivalent 
non-integrated-by-parts formulation: 


Find E V h such that for all v/ h E V h 

(3.10) B(v\w h ) GA L = 0 


where 

B(y, w)gal = / f w-(u t + V x .(y)) dxdt 
Ji n Jn 

+ / w(t") • [u]|+ dx 

J Q 

+ f X! f \ Mxt •h d (v(x_),v(x + );n)dxd< 

+ f f (( w ))x! [ f ( v ; n )]*! dxdt 


where h d denotes the flux dissipation term incorporated into the total numerical flux. 

Theorem 3.4. DG Global Entropy Norm Stability (Nonlinear Hyperbolic System). The 
variational formulation (3.10) for nonlinear systems of conservation laws with convex entropy extension and 
symmetric mean-value flux dissipation 


h SMv( v -.v+;n) = |A|smv[v]*+ , \A\ S mv = f |i(v(0);n)| j4 d0 

Jo 


is entropy norm stable with the following global balance: 
at- 1 


(3.11) 


(III Mil 1114 + + /"(<">* = (m)<k 

n= 0 \ e ££ / 


|I(n)| = [ 2(1-9) (i + (v(0); n) io - i'(v(l - 0); n)^) 


dd 


with 
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Proof. Repeated here from [2, 1]. Construct the energy balance for the interval — U n=o I n by 

setting w = v and evaluating the various integrals. Consider the time derivative integral 



v T u t dtdx 



U ti dt dx 




dx 


and combine with the jump integral across time slabs. From Corollary 3.2 

JJ v T u t dtdx + I v^Mufi dx = J n Mir dr + i||| Mg IH^ • 

When summed over all time slabs, the first term on the right-hand-side of this equation vanishes except for 
initial and final time slab contributions. Next, consider the spatial operator term and apply the divergence 
theorem 

(3.12) [ [ v T r Xi dxdt= f [ F' x . dx dt = [ [ - [F(v; n)J*+ dx dt 

7/"7n ’ J i n J n ■' l "e€£'^ e 


where F(v; n) = n, F'(v). Combining all the space terms and applying Corollary 3.3 

/4-pace S lx l (- [F(v; n)) l: + ((V»*+ ■ [f(n))l: + \ [v£+ • h“) dx dt 

= lx l \ M*! • ( hd + jf V ~ Mx! fcdt ■ 

In summary, collecting all terms and summing over time slabs we have 

N - 1 

S(v,v)gal = X! 

n=0 

= [ U(t?)dx- [ U(t°_)dx+'f i ( 1 - |||[v];l + ||| 2 ^ n + //r p ace) ■ 

Jn J n n=0 / 

When written in this form, it becomes clear that a sufficient condition for energy stability is that for all time 
intervals J" 



(3.13) 


//space > 0 


which serves as a design condition for the flux dissipation. 

//space = IT l \ M*- ' ( hd + l! {1 ~ W) [V] *- ^ 

= lx l \ M*t • ( hd + [ a - 9) Aitrm m:: <») *= * 

dedxdt 


The choice 
(3.14) 


= hg MV = [ |i(v(fl); n)| io [v]*+ dB 
J 0 


II n 

x '‘space 


2 £ ( (1 


9) A + (v(ey,n) Ao - 8 A -(v(ey,n) Ao ) d6 [v] x x + _ dxdt 


yields 
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= / y;/i[v]^-2 f\l-6) (A + (noy,n) Ao -A-{ni-0);n) Ao )de[w]l + _dxdt 

J/- eG(f 7e 2 Jo 

= \ lA“)l[vC! dxdt >0 . 

•' l ’'e€£ Je 

This completes the sufficient condition for energy decay in time. ■ 

An important observation to be made concerning (3.13) and energy stability is that any flux dissipation h d 
for which 


(3.15) 


[v]t * hg MV < [v]t * h d 


is also energy stable with increased energy decay. This observation is used in later sections in designing 
simplified flux functions for the DG method. 

3.2. DG Nonlinear Energy Analysis for Kinetic Boltzmann Moment Closure Hierarchies. 
In this section, we analyze nonlinear energy properties of the discontinuous Galerkin method assuming the 
kinetic Boltzmann moment closure structure discussed in Sect. 2.1. Recall that in this framework 

(3.16) u — (m exp(v • m)) , f l = {v, m exp(v • m)> 
with conjugate entropy given by 

(3.17) W(v) = (exp(v-m)) . 

Using these definitions, we briefy examine energy stability of the DG method for these moment closure 
hierarchies. 

Theorem 3.5. DG Global Entropy Norm Stability of Kinetic Boltzmann Moment Closure 
Hierarchies. The variational formulation (3.10) for nonlinear systems of conservation laws with convex 
entropy extension and kinetic Boltzmann moment closure symmetric mean-value flux dissipation 

hi SMV (v-,v + ;n) = |i| KS Mv[v]^ + , |A|ksmv = f (|t» • n| m<8> m exp(v(0) ■ m)) dO 

Jo 

is entropy norm stable with the following global balance: 

(3.18) 


; E (him!? ii6.,d + BM£>; S„ U{t - )dx = L U{t °- )dx 

n= 0 \ eeS / 


with 


|i(n)|= f 2(1 -9) (A + (v(0y,n) Ao - A~{v(l - Oy,n) Ao ) 


dO . 


Proof. By making the following generalizations 

C/(u) = ((v(u) • m) exp(v(u) - m) - exp(v(u) • m)) 
A 0 (v) = (m ® m exp(v • m)) 

A(v; n) = ((u • n) m 0 m exp(v • m)) 

A ± (v; n)^ o = ((u • n) 1 * 1 m ® m exp(v • m)) , 
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we can appeal once again to the space-time slab jump identities stated in lemma 3.1 and corollaries 3.7 
and 3.3. Using these definitions and results, the proof of theorem 3.4 applies without alteration up to and 
including the equation 


"s^ac, = /E / \ Mxt ■ ( hd + / l(1 - 20) A ' { ' m M*- dd ) dX dt 
= J ]E f\ M»! • (h d + J (1 - 6) ii(v(0)) [v]*+ de^j dxdt 

- 1 E l \ Mx! • [ * ^(vw) m:: dBdxdt 


and the design condition 
(3.19) 

In the present case, the choice 


''space > 0 - 


(3.20) 
is sufficient 


h‘ , = hf < 


KSMV 


= [ (|u • n| 
Jo 


m ® m exp(v(0) • in)) dd [v] 


*+ 

X- 


"space = I J2 1 \ Mx! • 2 l ((! - 8) ^ + (v(0); n)^ -flA- (v(fl);n) io ) d6 [v]** dxdt 

= dd[v}l+ dxdt 

= l E l \ w:: • ti(“)i Mx- dxdt>Q . 


e££ ' 

This completes the sufficient condition for energy decay in time for kinetic Boltzmann moment closure 
hierarchies. ■ 

4. Simplified Numerical Flux Formulas for the DG Method. The theoretical results of Sect. 3 
provide the framework for constructing, analyzing, and proving energy stability for a number of simplified 
numerical flux functions. This task is undertaken in the remainder of this section. We are unaware of any 
previous DG analysis for systems (m > 1) of nonlinear conservation laws which rigorously establishes energy 
stability for the fluxes considered here. Throughout, we use the notation A = R A R T as defined earlier with 
A = diag(Ai, . . . , A m ) assuming ordered entries Aj < . . . < A m . For numerical fluxes such as the SHHLE and 
SHLLEM flux, we also require that Ai and A m be distinct in order that the construction be well defined. 

• Symmetric Lax- Friedrichs Flux (SLF) 


(4.1) 

with 


h S LF(v-,v+;n) 


\ (f(v_;n) + f(v + ;n)) - [u(v)]^ 


A max — SUp (A m (v(^))) . 

0 < 9<1 

• Symmetric Lax-Friedrichs Matrix Flux (SLFM) 

(4.2) hsLFM ( v -) v +;n) = - (f(v_;n) +f(v + ;n)) - - hs LFM (v_, v+; n) 

with 

^SLFm( v -> v +j n ) = ^max ^ (^o( V ~) + ^o( v +)) [ y ] X - 
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and 


^max — Slip (A m (v(0))) 

O<0< 1 


• Symmetric Harten-Lax-van Leer Flux (SHLLE) 

(4.3) hsHLLE(v-, v+;n) = - (f(v_;n) 4* f(v+;n)) - - hs HLLE (v_, v + ; n) 


with 


+ A n 


2A n 


, d < \ ^tnax "T rr/ \i-{- ^^max^min + 

^SHLLE ( v — y v + » n ) = T — T [ f ( v ; n )J- " T L U ( V )]- 

Amax '''min '''max / 'min 


and 


A 


max 


sup max(0, A m (v(0))) , A min = inf min(0, Ai(v(0))) . 

0<9<1 0<^<1 


• Symmetric Harten-Lax-van Leer Modified Flux (SHLLEM) 

(4.4) hsHLLEM(v_, v + ;n) = - (f(v_;n) + f(v + ;n)) - - hg HLLEM (v_, v+; n ) 

with 


/ \ A max *+• A m j n . . . .1 2A max A m i n r . i f . . _ - r-i-f- 

(v— , v+; n) = [f (v; n)]I - - — [u(v)]+ - ^ / fi{6) n ® n [v]I 

^max '''min ^max ''min ^_2 JO 


h? 


SHLLEM \ 


dd 


with 

max(0,A m (fl)) +min(0,Ai(fl)) _ 2 max(0,A m (fl)) min(0,Ai(fl)) _ . 

max(0, A m (0)) - min(0, Ai (0)) * max(0, A m (0)) - min(0, Ai(0)) 

and 

Amax - sup max(O,A m (v(0))) , A min = inf min(0, Ai(v(fl))) . 

O<0<1 0<*<1 


• Discrete Kinetic Symmetric Mean- Value Flux (DKSMV) 


(4-5) hoKSMV (v_ , v+ ; n) = ^ (f(v_;n) + f(v + ;n)) - i h& KSMV (v_,v+;n) 

with 

h DKSMv( v -> v +; n ) = \ (|v n|m® m(exp(v_ • m(u)) + exp(v+ • m(u))} [v]± 

= ^ (|v • n\ m 0 m exp(v_ • m(u))) [v]i 
-f ^ (\v • n| m 0 m exp(v+ * m(u))) [v]l 

Observe that explicit path ^-integration has been avoided in all these simplified fluxes (except correction 
terms in SHHLEM). In addition, we have the following theorem: 

THEOREM 4.1. Energy Stability of Simplified Flux Formulas The variational formulation (3.10) 
for nonlinear systems of conservation laws utilizing any of the candidate approximate numerical fluxes (4-1), 
(4-2), (4-3), (4-4) f (4-5) * s entropy norm stable in the sense of Theorem 3.4 or 3.5. 

Proofs: Given on a cases-by-case basis. 
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Symmetric Lax-Friedrichs Flux (SLF): 


Mi! -hsMv( v -i v +; n ) = Mi! • [ l^(v(0);n)| ii [v]*+<» 

Jo 

= [ Mi! • A(y(0); n) |A(v(0))| R T (y(ey, n) [v]*+ d6 
Jo 

< sup (A m{v( 6 ))) [ Mi! •fl(v(«);n)f(v(S)!n)[v]'+ dO 
O<0<1 Jo 

= sup (A m (V(0))) f Mi! • A>(v(0)) M*! dd 
O< 0<1 Jo 

= sup (A m (v(0)))[v]*+ -[u(v)]*+ 

O<0<1 

= Mi! - h sLF( v -» v +;n) ■ 


Symmetric Lax-Friedrichs Matrix Flux (SLFM): 


M^! -^Mv(v-,v + ;n) = [v]£ • [' |A(v(0); n)\ Ao [v]»+ dS 

Jo 

= f Mil • A(v(0); n) |A(v(0))| R T (V(0); n) [v]£ d6 
Jo 

< sup (A m (v(0))) f Mi! • R(v(e);n)R r (V{ey i n) [v]* + d 6 
O<0<1 Jo 

= sup (A m (v(0))) f M! + • Mr dS 

0 < 9 <l Jo 

Examining the scalar function 

9 ( 6 ) = [v]± • A o (v(0))[v]t = [v]i • ^(v(0))[v]t , 

differentiation yields 

5" = 0(vW)[M!,Mt,[v]!,[v]+] 

where the right-hand-side term denotes the rank-4 contraction to a scalar. But for systems derived as 
moments of the kinetic Boltmann equation, we have from Sect. 2.1 that 


d A U 

dv 4 


[z, z, z, z] > 0, |z| / 0 . 


Consequently, g( 6 ) is convex for all 0 , 


9(0)<(l-6)9(0)+0g(l) , 

so that 

Mi! • [ m:! • 4>(v(*)) Mi! m<\ (M£ • Miy-) M*! + Mil • 4>(v+) Mil) 

= \ Mi! • (^(v_) + io(v+)) Mi! ■ 

This yields 

Mi! ' hsMv( v -> v+; n) < sup (A m (v(0))) \ [v]^ • (a 0 (v_) + i 0 (v + )) [v]£ 
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“ Mxl ‘ ^SLFM ( v - i V +5 n ) 


Symmetric Harten-Lax-van Leer Fluxes (SHLLE) and (SHLLEM): 

Our first task will be to prove that 

M± * hsMv( v -> v +5 n ) < M- ' ^SHLLEM ( v --^ v + 5 n ) 

followed by 

Ml ’ ^shLLEm( v -i v +i n ) ^ Ml * ^SHLLe( V --i v +5 n ) • 

To do so, consider the symmetrization matrices A and .4 q written in the form (2.37) 


A = Xj fj ® fj Ap = ^ fj ® fj 

<=i i=i 

find the following useful 2x2 matrix identity for A £ 1R 2x2 and Ao € 1R 2x2 
(4.6) 


, _ max(0, A 2 ) + min(0, Xy) - 2 max(0, A 2 ) min(0,Ai) - 

— 77, x \ : 77, x 7 A ' 77T~7 \ /n \ \ ^0 


max(0, A 2 ) - min(0, Ai) max(0, A 2 ) — min(0, Ai) 

as can be easily verified by substitution. Generalizing to m x m matrices with Ai < A 2 < . . . < A m and 
Ai < A m , we can only represent the extremal values of A* exactly using the (4.6) ansatz whenever Ai A m < 0. 
Consequently, we have a slightly more complicated identity for general m > 2 

, 7 , _ max(0, A m ) + min(0, Aj ) r 2 max(0,A m ) min(0 , AQ - _ 

'Ao mayfO — min/Tl mayffl A „ "l — minff) ° ' f * * * 


with 


fi = 


max(0, A m ) — min(0, Ai) max(0, A m ) — min(0, Ai) 

max(0, A m ) + min(0, Ai) 2 max(0, A m ) min(0, Ai) 
■A* — 


t=2 


-N 


max(0, A m ) — min(0, Ai) * max(0, A m ) - min(0, Ai) 

_ max(0, A m ) (A~ - min(0, AQ) + min(0, A t ) (At - max(0, A m )) > 
max(0. A m ) — min(0, Ai) “ 

Next, consider the local path integral form of hg MV 

h SMv( v -> v +; n ) = [ l^(v(0)U o [v]td0 

Jo 

( m -i \ 

<n (*)£(?(*); n) - a 2 (9)A 0 (v(9)) - £ /i(0)f, <8> r< [v]t d9 


=/ 


i=2 


with 


and 


rnax(0,A m (fl)) + min(0, Ai(fl)) (g , _ 2 max(0,A m (g)) min(0,Ai) _ 

* max(0, X m {9)) — min(0, \i(9)) ' max(0, A m (0)) — min(0, Aj(0)) * ~ 

max(0, A m (v(0))) + min(0, Ai(v(0))) 


cn{9) = 
a 2 (9) = 


max(0, A m (v(0))) — min(0, Ai(v(0))) 
2max(0, A m (v(fl))) min(0, Ai(v(fl))) 
max(0, A m (v(0))) - min(0, Ai(v(0))) 
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In addition, we will define the perturbed ratios 2 

- m _ (max(0, A ro (v(fl))) + 6(0)) + (min(0, A t (v(g))) - 7(g)) 

<TlW (max(O,A m (v(0))) + 5(d))-(min(O,A 1 (v(0)))- 7 W) ’ 

- (f) s _ 2(max(0, A m (v(fl))) + 5(0)) (min(0, A t (v(fl))) - jrW) 
a2( ’ (max(0,A m (V(*))) + *(*)) - (min(0, \i(V(9))) - 7 (*)) 

for nonnegative bounded functions 5(9) > 0 and 7 (0) >0,06 [0, 1]. Examination of the scalar quantity 

II = J [v]i • (<ri(0)i4(v(0);n) - 5 2 (0)i o (v(0))) Mi dO 
- J\v]± • (a 1 (9)A(V(9)-, n) - <r 2 (0)i o (v(0))) [v]t d9 
= J\y ]tR(9) • [(MO) - MO))HO) - (MO) - MOWmxm) R T (0) Mi d9 

reveals that II >0 since for each component A * of A, i = 1 , . . . , m, (omitting the dependence on 9) 

r _ \ \ _ (~ _ \ _ 2 7 max(Q, A m ) (max(0, A m ) + 6 - \j) 

o\ °\ * o 2 o 2 (( max (0 ,A m ) + 5) - (min(0,Ai) -7)) (max(0,A m ) - min(0,Ai)) 

2 <5 min(0, Ai) (min(0, Ai) — 7 — Aj) 

((max(0, A m ) + 5) - (min(0, Ai ) - 7)) (max(0, A m ) - min(0, Ai )) 

> 0 . 

Define infimum and supremum values of min(0, Ai (0)) and max(0, A m (0)) respectively in the interval 0 E [0, 1] 
as 


A max = sup max(0, A m (v(0))) , A min = inf min(0, Ai (v(0))) 
O<0<1 o<»<i 

and set c(9) and (5(0) as follows 

e(0) = A max - max(0, A m (v(0))) > 0 , <5(0) = min(0, Ai(v(0))) - A min > 0 

This renders 07 (0) and a 2 (9) now 0-independent, i.e. 

— ^max "1“ A m j n — 2A max A m i n 
cq = t ZTT~ ’ a2 

^max / 'min 

Consequently, for II > 0 and /*(0) > 0 we have 


^max A m i n 


A- ■ h SMV = M- • f (cri(0)A(v(0);n) - ct 2 (0)A o (v(0))) [v]1 d9 

r\ m -1 

-Mi • / /<(*) Mi de 

J ° t=2 

< Mi • J (<?i(0)A(v(0);n) - 5 2 (0)i o (v(0))) [v]i d9 

rl m -1 

-Mi • / 5Z f^o)fi ® fi Mi d o 

*'° i~2 

= Mi • ( 5 i J ^(v(0); n) Mi do -a 2 J i o (v(0)) Mi do) 


2 Our strategy will be to later define the nonnegative function 6(9) as the “gap” between local value of max(G, A m (#)) and 
the supremum value in the interval $ 6 [0, 1] and similarly 7(0) will represent the gap between the min(0, Ai(0)) and the infimum 
value. 



16 


m-l 


-M- • [ /<(*) fi ® fi M- dd 

Jo i= 2 

= [v]t ■ (ax [f(v;n)]t - [v]t • 5 a [u(v)]i) 

-1 m-l 

-Ml • / V /,(6>) fi ® fi Ml d0 

Jo i= 2 

= Ml ' hsHLLEM 

< Ml ' ( 5 i [f(v;n)]l - [v]l • a 2 [u(v)]l) 


— Ml ' hsHLLE 


which completes the proof for the SHLLE and SHLLEM fluxes. 
Discrete Kinetic Symmetric Mean- Value Flux (DKSMV): 


Ml ■ h KSMV = Ml ■ f (|u • n| m ® m exp (v(0) • m(v))) Ml dO 

= [v]± • • n\ m ® m J exp (v(0) • m(u)) dO ^ [y]t 

= (^\v • n| ([v]t • m) 2 J exp (v(0) • m(u)) dO ^ . 

Considering the scalar function 

9(0) = exp (y(9) • m(v)) 

followed by twice differentiation 

9"(0) = ([v]i • m(t?)) 2 exp (v(0) • m(f)) > 0 . 

Hence, g{9) is yet smother convex function so that 

9(0) <(1-0) 9(0) + 0g(l) 


and finally 


M- 


a KSMV 


— r v i+ 


• f (|u * n| m ® m exp (v(0) * m(u)) ) [v]l dO 

Jo 


< [v]t • -(|u ■ n| m 0 m (exp(v_ • m(u)) + exp(v+ • m(u)))) [v]t 
= [vjl * hi} KSMV ■ 


5. Concluding Remarks. The analysis of this paper confirms energy stability for several numerical 
flux functions that are of practical merit when used in computational fluid dynamics computations. Even 
so, the theoretical framework developed here applies more generally and has application to many nonlinear 
conservation law systems with entropy extensions that are not explicitly discussed here. In these settings, 
the analysis presented may be invaluable because there may not be the large body of numerical methods 
developed before hand to guide the development of new numerical fluxes, discretization, and stabilization. 
Our general goal is to pursue these new problem areas in forthcoming work. 


Appendix A. The exp(v m(v } I)) Boltzmann Moment Structure for a 7 -Law Equation of 
State. For a 7-law (polytropic) gas, one has 


P = (7 - 1 )P«» T={ 7 - l)f • 
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Following Perthame [18], we consider the following Maxwellian in JR d for a 7 -law gas: 


(A.l) 

with 

and 


f{p,u,T\vJ) = 


a(7,d) 


e -(\u~ v \*/2+I*)/T 


8 = 


__1 d 

7-1 2 


a( 7 ,d) = f e ^ 2 dv • f e 
JlR d J IR + 


dl 


Using this particular form, Perthame shows that the Euler equations for a 7 -law gas are obtained as the 
following moments 


(A. 2 ) 


1 

m (v,I) = | v 

>| 2 / 2 + I s 


u=(m /), r = (i>iin/), {■) = f [ (-) dl dv . 

J IR d J 1R+ 

The nonobvious energy moment |u| 2 /2 + I s was devised by Perthame rather than the more standard moment 
|t>| 2 /2 + 1 (see for example [6, 7]) in order that a classical Boltzmann entropy H(f) = / log / be obtained. 

Let us now verify that this choice of moments yields an exponential form for the conjugate entropy 
function of the form 

^( v ) = (/) = (c(7,d)exp(v • m(u,/))) , c(7,d)>0 

which is sufficient for our purposes. Inserting the expression for 6 into the temperature term appearing in 
the Maxwellian yields 

1 P 


(A.3) 


f(p,u,T-,v,I) = 


-(|«-«| a /2 +I l )/T 


a(j,d) TV( 7 - 1 ) 

and compare this with the expression exp(v • m) obtained using the entropy function 3 


U( u) = - 


ps 


so that 


2_ 4- _2 bill ’ 

_ ttt _ I I - 1 + r 1 2T 

Ji 

T 


v = = 


(7-1) 

log ~ 


T 1 

T 


and finally 


exp(v • m(v, /)) = e 


- p"7/(7-l) 


Q -(\u-v\ 2 /2+I*)/T 


yi/(7-l) 

Comparing with (A.l), we obtain the exponential form for the conjugate entropy function 

W(v) = (/) = (c( 7 ,d) exp(v • m(u, /))) 

with 

e 7/(7-l) 

c(7, d) = - > 0 . 

a(l,d) 


3 The choice of 1/(7 — 1) scaling of the entropy function comes from our desire to match the p/T 1 /^-*) term appearing in 
Perthame’s Maxwellian 



